data_2 %<>%
  mutate(income_recode2 = case_when(q26_survey_income_recode %in% c(1, 2) ~ 1,
                                    q26_survey_income_recode == 3 ~ 2,
                                    q26_survey_income_recode == 4 ~ 3,
                                    q26_survey_income_recode == 5 ~ 4,
                                    q26_survey_income_recode == 6 ~ 5,
                                    q26_survey_income_recode %in% c(7, 8) ~ 6,
                                    TRUE ~ NA_real_),
         q22_survey_working_recode = case_when(q22_survey_working == "Ya" ~ 1,
                                               q22_survey_working == "Tidak" ~ 0,
                                               TRUE ~ NA_real_),
         q23_survey_look_work_recode = case_when(q23_survey_look_work == "Ya" ~ 1,
                                                 q23_survey_look_work == "Tidak" ~ 0,
                                                 TRUE ~ NA_real_)
  )

test_mod1a <- lm(income_recode2 ~ pass_final, data = data_2 %>% filter(abs(forcing_final) < 1))
test_mod2a <- lm(q22_survey_working_recode ~ pass_final, data = data_2 %>% filter(abs(forcing_final) < 1))
test_mod3a <- lm(q27_survey_job_satis ~ pass_final, data = data_2 %>% filter(abs(forcing_final) < 1))

observations <- c(nobs(test_mod1a), nobs(test_mod2a), nobs(test_mod3a))

test_mod1a <- coeftest(test_mod1a, vcov=vcovHC(test_mod1a,type="HC0"))
test_mod2a <- coeftest(test_mod2a, vcov=vcovHC(test_mod2a,type="HC0"))
test_mod3a <- coeftest(test_mod3a, vcov=vcovHC(test_mod3a,type="HC0"))

#Javanese advantage table

table <- list(test_mod1a, test_mod2a, test_mod3a)

note_text <- paste(" Beta coefficients from OLS regression. Standard errors were calculated using the Huber-White (HC0) correction. 
                   The outcomes measure are indexed values capturing (1) Javanese preferentialism among Javans and (2) among non-Javans, (3) regional preferentialism,
                   (4) religious resentment, (5) perceptions of corruption, (6) national identification.")

table = stargazer(table, 
                  type = 'latex', 
                  title = "The Effect of Passing Specialist Competence Examination (SKB) On Work Outcomes",
                  label = 'tab:service_effect_work',
                  model.names = F,
                  model.numbers = T,
                  digits = 2,
                  column.separate = c(1, 1, 1),
                  column.labels = c("Income (m, IDR)", "Employed", "Job Satisfaction (1-4)"),
                  multicolumn = T,
                  dep.var.labels = NULL, 
                  add.lines = list(#c("Subset", "----", "Non-CS", "---", "Non-CS", "---", "Non-CS"),
                    c("Observations", observations),
                    c('Bandwidth', rep(c('1\\%'), 3))),
                  covariate.labels = c("Passed Test"),
                  star.cutoffs = c(0.05, 0.01, 0.001),
                  #float.env = 'sidewaystable',
                  keep.stat = c("n"),
                  notes = NULL,
                  notes.align = 'l')


write_latex(table[-c(10, 11, 12, 18, 21, 27)], note_text, './_4_outputs/tables/table_a16.tex', .8)

